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Abstract 

We review results obtained using Shell Model Monte Carlo 
(SMMC) techniques. These methods reduce the imaginary-time 
many-body evolution operator to a coherent superposition of one- 
body evolutions in fluctuating one-body fields; the resultant path 
integral is evaluated stochastically. After a brief review of the 
methods, we discuss a variety of nuclear physics applications. 
These include studies of the ground-state properties of p/-shell 
nuclei, Gamow- Teller strength distributions, thermal and rota- 
tional pairing properties of nuclei near N = Z, 7-soft nuclei, and 
/3/3-decay in ^''Ge. Several other illustrative calculations are also 
reviewed. Finally, we discuss prospects for further progress in 
SMMC and related calculations. 
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1 Introduction 

The notion of independent particles moving in a common one-body po- 
tential is central to our description of atoms, metals, and hadrons. It is 
also realized in nuclei, and the shell structure associated with the magic 
numbers was first put on a firm basis in 1949, when the magic num- 
bers were explained by an harmonic oscillator spectrum with a strong, 
inverted (with respect to the atomic case) spin-orbit potential [Q. 

But nuclei differ from the other quantal systems cited above in that 
the residual interaction between the valence fermions is strong and so 
severely perturbs the naive single-particle picture. This interaction 
mixes together many different configurations to produce the true eigen- 
states and, because of its coherence, there emerge phenomena such as 
pairing, modification of sum rules, deformation, and collective rotations 
and vibrations. An accurate treatment of the residual interaction is 
therefore essential to properly describe nuclei. 

A hierachy of models has been developed to treat various aspects of 
the mean field and correlations; Skyrme Hartree-FockQ for the mean 
field, Hartree-Fock Bogolioubov for the pairing correlations , and RPA 
B for particle-hole correlations. But short of a complete solution of the 
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many-body problem, the shell model is regarded as the most fundamen- 
tal approach for studying nuclear structure, provided that the single- 
particle space in which calculations are performed is large enough. Suc- 
cess was demonstrated early in the Op-shell [|], in the Os-ld in the 
0/7/2 region Q, and more recently in the lpO/-shell using interactions 
developed by Kuo and Brown Q and modified by Zuker and Poves 

The traditional numerical approach to the nuclear shell model is to 
diagonalize the hamiltonian matrix within a limited many-body basis 
]lOt . Direct diagonalization gives the wave functions which can then 
be used to calculate the properties of specific levels. The shell model 
treated in this way has been quite successful in describing sd-shell nuclei 
and has also been used to successfully describe various properties 
of the lower p/-shell ||l2|, However, the direct-diagonalization pro- 
cedure is limited by memory requirements that grow combinatorially 
with the size of the single-particle space, Ns, and with the number of 
valance particles, Ny In the pf -shell where ^°Ca is the core, the largest 
systematic study thus far treats A = 48,49 an achievement that 
required approximately five generations of technological and software 
improvement. 

To circumvent the limitations of direct diagonalization, we have pro- 
posed an alternative treatment of the shell model 0, ^ [l^ . 
Our methods are based on a path-integral formulation of the imaginary- 
time many-body propagator, exp(— where H is the shell model 
hamiltonian, and /3 is the reciprocal of the temperature. We employ 
the Hubbard-Stratonovich (HS) transformation to recast the two- 
body terms in the exponential into one-body terms with an integration 
over auxiliary fields. Thus the nucleons are seen as independent parti- 
cles moving about in fluctuating fields Such auxiliary field meth- 
ods have also been applied to condensed matter systems such as the 
Hubbard model |2^, , yielding important information about electron 
correlations and magnetic properties. 

Having circumvented the diagonalization of the hamiltonian, we are 
left with a large multidimensional integration to perform, for which we 
use the Monte Carlo sampling techniques of Metropolis et al . Real- 
istic nuclear hamiltonians often have a sign problem that hampers the 
Monte Carlo quadrature. We have overcome this by an extrapolation 
procedure from a family of hamiltonians, all of which have no sign prob- 
lem and are close to the full hamiltonian. It is important to note that 
these calculations are carried out on large scale parallel computers, a 
technology that is essential to implement these methods. 

In this article, we briefly outline SMMC methods in Section II, then 
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concentrate our discussion on the physics resuhs obtained. We discuss 
in Section III our major results, which include ground state properties of 
p/-shell nuclei, Gamow- Teller strengths and distributions, thermal and 
rotational pairing properties of nuclei near N = Z, 7-soft nuclei, and 
/3/3-decay. In Section IV, several other illustrative studies are reviewed. 
We conclude in Section V with a discussion of future prospects. Our 
coverage of the field extends through February 1997. 



The nuclear shell model is defined by a set of spin-orbit coupled single- 
particle states with quantum numbers Ijm denoting the orbital angular 
momentum (l) and the total angular momenta (j) and its z-component, 
m. Although non-spherical one-body potentials are a common efficiency 
used in describing deformed nuclei, for the rotationally invariant hamil- 
tonians used in SMMC so far, these states have energies eij that are 
independent of m. The single-particle states and energies may be differ- 
ent for neutrons and protons, in which case it is convenient to include 
also the isospin component ^3 = ±1/2 in the state description. We 
will use the label a for the set of quantum numbers Ijm or Ijmt^, as 
appropriate. 

In the following we briefly outline the formalism of the SMMC method. 
We begin with a brief description of statistical mechanics techniques 
used in our approach, then discuss the Hubbard-Stratonovich transfor- 
mation, and end with a discussion of Monte Carlo sampling procedures. 
We refer the reader to previous works jlj, |l^ for a more detailed expo- 
sition. 

2.1 Observables 

SMMC methods rely on an ability to calculate the imaginary-time many- 
body evolution operator, exp(— where /3 is a real c-number. The 
many-body hamiltonian can be written schematically as 



where O is a density operator, V is the strength of the two-body in- 
teraction, and £ a single-particle energy. In the full problem, there are 
many such quantities with various orbital indices that are summed over, 
but we omit them here for the sake of clarity. 



2 Methods 




(1) 
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While the SMMC technique does not result in a complete solution to 
the many-body problem in the sense of giving all eigenvalues and eigen- 
states of H, it can result in much useful information. For example, the 
expectation value of some observable fl can be obtained by calculating 

Here, f3 = is interpreted as the inverse of the temperature T, and 
the many-body trace is defined as 

TrX^^(z|X|z) , (3) 

i 

where the sum is over many-body states of the system. In the canonical 
ensemble, this sum is over all states with a specified number of nucle- 



ons (implemented by "number projection" |I5, |l8j), while the grand 
canonical ensemble introduces a chemical potential and sums over all 
many-body states. 

In the limit of low temperature (T — > or /? — > oo), the canonical 
trace reduces to a ground state expectation value. Alternatively, if |<I>) 
is a many-body trial state not orthogonal to the exact ground state, j^*), 
then e~^^ can be used as a filter to refine |$) to j^") as (3 becomes large. 
An observable can be calculated in this "zero temperature" method as 



(4) 



If Q is the hamiltonian, then (Q) at /3 = is the variational estimate 
of the energy, and improves as /? increases. Of course, the efficiency of 
the refinement for any observable depends upon the degree to which |<I>) 
approximates \^>). 

Beyond such static properties, e"^^ allows us to obtain some infor- 
mation about the dynamical response of the system. For an operator 
f2, the response function RniT) in the canonical ensemble is defined as 

Rn{r) ^ ^^Jl; ^ {n\r)m), (5) 

where ^^t) = e^^VL^ e^'^^ is the imaginary-time Heisenberg operator. 
Interesting choices for 17 are the annihiliation operators for particular 
orbitals, the Gamow- Teller, Ml, or quadrupole moment, etc. Inserting 
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complete sets of A-body eigenstates of H |/)} with energies -Ej,/) 
shows that 

Rn{r) = ^Y.^-''''\mi)?e-^^'''-'''\ (6) 

if 

where Z = X^^e"''^* is the partition function. Thus, Rn{T) is the 
Laplace transform of the strength function Sn{E): 

/oo 
e-^''Sn{E)dE ; (7) 
-OO 

Sn{E) = |^e-'^^'|(/|^lH)r5(i;-iS/ + £;,). (8) 

ft 

Hence, if wc can calculate Rq(t), Sii{E) can be determined. Short of 
a full inversion of the Laplace transform (which is often numerically 
difficult), the behavior of i?n(r) for small t gives information about the 
energy- weighted moments of Sn- In particular, 

/oo 1 
Sa{E)dE = -Y.e-^^'^lifm^ = {n^^)A (9) 
-oo ^ i 

is the total strength, 

/oo 1 
S^{E)EdE = J2^-f^'''\ifm'iEf - E,){10) 
^ if 

is the first moment (the prime denotes differentiation with respect to 

It is important to note that we usually cannot obtain detailed spec- 
troscopic information from SMMC calculations. Rather, we can calcu- 
late expectation values of operators in the therniodynaniic ensembles 
or the ground state. Occasionally, these can indirectly furnish prop- 
erties of excited states. For example, if there is a collective 2+ state 
absorbing most of the E2 strength, then the centroid of the quadrupole 
response function will be a good estimate of its energy. But, in general, 
we are without the numerous specific excitation energies and wavefunc- 
tions that characterize a direct diagonalization. This is both a blessing 
and a curse. The former is that for the very large model spaces of in- 
terest, there is no way in which we can deal explicitly with all of the 
wavefunctions and excitation energies. Indeed, we often don't need to, 
as experiments only measure average nuclear properties at a given exci- 
tation energy. The curse is that comparison with detailed properties of 
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specific levels is difficult. In this sense, the SMMC method is comple- 
mentary to direct diagonalization for modest model spaces, but is the 
only method for treating very large problems. 



2.2 The Hubbard-Stratonovich transformation 

It remains to describe the Hubbard-Stratanovich "trick" by which e^^^ 
is managed. In broad terms, the difficult many-body evolution is re- 
placed by a superposition of an infinity of tractable one-body evolutions, 
each in a different external field, a. Integration over the external fields 
then reduces the many-body problem to quadrature. 

To illustrate the approach, let us assume that only one operator O 
appears in the hamiltonian (|l]). Then all of the difficulty arises from 
the two-body interaction, that term in H quadratic in C If iJ were 
solely linear in O, we would have a one-body quantum system, which is 
readily dealt with. To linearize the evolution, we employ the Gaussian 
identity 



2tt 



/OO 
dae-i'^l^l^'e-'^''; h = eO + sVaO . (11) 
-OO 



Here, /i is a one-body operator associated with a c-number field ct, and 
the many-body evolution is obtained by integrating the one-body evo- 
lution Ua = e~^^ over all a with a Gaussian weight. The phase, s, is 
1 if < or i if F > 0. Equation ( pT| ) is easily verified by complet- 
ing the square in the exponent of the integrand; since we have assumed 
that there is only a single operator O, there is no need to worry about 
non-commutation. 

For a realistic hamiltonian, there will be many non-commuting den- 
sity operators Oa present, but we can always reduce the two-body term 
to diagonal form. Thus for a general two-body interaction in a general 
time-reversal invariant form, we write 

a a 

where Oa is the time reverse of Oa- Since, in general, [Oa,Op] ^ 0, we 
must split the interval f3 into Nt "time slices" of length A/3 = P/Nt, 

^-PH ^^^-ApH^m^ (13) 

and for each time slice n = l,...,Nt perform a linearization similar 



to Eq. 11 using auxiliary fields aan- Note that because the various 
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Oct need not commute, the representation of e^^^^ must be accurate 
through order to achieve an overall accuracy of order A/3. 

We are now able to write expressions for observables as the ratio of 
two field integrals. Thus expectations of observables can be written as 



where 



^ n—1 a ^ 



and 



h„^^ (e* + SaVaaan) O a + (Sa + SaVaCr*„) Oa ■ (16) 

This is, of course, a discrete version of a path integral over a. Because 
there is a field variable for each operator at each time slice, the dimension 
of the integrals T>a can be very large, often exceeding 10^. The errors 
in Eq. |lj are of order A/3, so that high accuracy requires large Nt and 
perhaps extrapolation to TVt = oo (A/3 = 0). 

Thus, the many-body observable is the weighted average (weight 
Wa) of the observable fla calculated in an ensemble involving only the 
one-body evolution Ua- Similar expressions involving two a fields (one 
each for e^"^^ and e^^^^'^'^^) can be written down for the response func- 
tion (^), and all are readily adapted to the canonical or grand canonical 
ensembles or to the zero-temperature case. 

An expression of the form (|lj) has a number of attractive features. 
First, the problem has been reduced to quadrature — we need only cal- 
culate the ratio of two integrals. Second, all of the quantum mechanics 
(which appears in fl^) is of the one-body variety, which is simply han- 
dled by the algebra of x Ng matrices. The price to pay is that we 
must treat the one-body problem for all possible a fields. 



2.3 Monte Carlo quadrature and the sign problem 

The manipulations of the previous sections have reduced the shell model 
to quadrature. That is, thermodynamic expectation values are given 
as the ratio of two multidimensional integrals over the auxiliary fields. 
The dimension D of these integrals is of order N'^Nt, which can exceed 
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10^ for the problems of interest. Monte Carlo methods are the only 
practical means of evaluating such integrals. In this section, we review 
those aspects of Monte Carlo quadrature relevant to the task at hand. 
We begin by recasting the ratio of integrals in Eq. (|lj) as 

(n) = J d^aPM. , (17) 

where 

bmce / d^aPa^ = 1 and P^ > 0, we can think of Pg. as a probabil- 
ity density and (il) as the average of fta weighted by P^- Thus, if 
{cts, s = 1, . . . , S"} are a set of S field configurations randomly chosen 
with probability density Pa, we can approximate (fl) as 

1 S 

, (19) 

where ils is the value of fi^ at the field configuration <7s- Since this 
estimate of (il) depends upon the randomly chosen field configurations, 
it too will be a random variable whose average value is the required 
integral. To quantify the uncertainty of this estimate, we consider each 
of the Cls as a random variable and invoke the central limit theorem to 
obtain 

4^) = ^ / d^'aPaiQa - « ^ piQs - . (20) 

This variance varies as S^^^^. 

We employ the Metropolis, Rosenbluth, Rosenbluth, Teller, and 
Teller algorithm to generate the field configurations Og, which 

requires only the ability to calculate the weight function for a given 
value of the integration variables. This method requires that the weight 
function Wa- must be real and non-negative. Unfortunately, many of 
the hamiltonians of physical interest suffer from a sign problem, in that 
Wa- is negative over significant fractions of the integration volume. To 
understand the implications of this, let us rewrite Eq. dlTl ) as 

(n) = / d^aPa-i-a^a , (21) 
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where 

\W,\ 

and — Wa/ I Wa \ is the sign of the real part of Wa. (Note that 
since the partition function is real, we can neglect the imaginary part.) 
Since | \ is non-negative by definition, we can interpret it, suitably 
normalized, as a probability density, so that upon rewriting (p7|) as 



we can think of the observable as a ratio in which the numerator and 
denominator can be separately evaluated by MC quadrature. Leaving 
aside the issue of correlations between estimates of these two quantities 
(they can always be evaluated using separate Metropolis walks), the 
fractional variance of (J7) will be 



in) - V + (<i>)2 ' ^^"^^ 

which becomes unacceptably large as the average sign ($) approaches 
zero. The average sign of the weight thus determines the feasibility of 
naive MC quadrature. In most cases ($) decreases exponentially with 
P or with the number of time slices |2J]. 



It has been shown |15 that for even-even and N ~ Z nuclei there 
is no sign problem for hamiltonians if all Vq < 0. Such forces include 
reasonable approximations to the realistic hamiltonian like pairing plus 
multipole interactions. However, for an arbitrary hamiltonian, we are 



not guaranteed that all Vq < (see for example |l^, 15 ). However, we 
may expect that a realistic hamiltonian will be dominated by terms like 
those of the schematic pairing plus multipole force (which is, after all, 
why the schematic forces were developed) so that it is, in some sense, 
close to a hamiltonian for which the MC is directly applicable. Thus, the 
"practical solution" to the sign problem presented in Ref. |jl^ is based 
on an extrapolation of observables calculated for a "nearby" family of 
hamiltonians whose integrands have a positive sign. Success depends 
crucially upon the degree of extrapolation required. Empirically, one 
finds that, for all of the many realistic interactions tested in the sd- and 
p/-shells, the extrapolation required is modest, amounting to a factor- 
of-two variation in the isovector monopole pairing strength. 
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Based on the above observation, it is possible to decompose H in 



The "good" hamiltonian Hq includes, in addition to the one-body terms, 
all the two-body interactions with Va < 0, while the "bad" hamiltonian 
Hb contains all interactions with Vq > 0. By construction, calculations 
with Hq alone have = 1 and are thus free of the sign problem. 

We define a family of hamiltonians Hg that depend on a continuous 
real parameter g as Hg = f{g)Ha + ff^^s, so that i?g=i = H , and f{g) 
is a function with /(I) — 1 and f{g < 0) > that can be chosen to 
make the extrapolations less severe. (In practical applications f{g) = 
1 — (1 — 5)/x with X « 4 has been found to be a good choice.) If the 
Va that are large in magnitude are "good," we expect that Hg^Q — Hq 
is a reasonable starting point for the calculation of an observable (il). 
One might then hope to calculate {n)g = Tr (f7e"'^^«)/Tr (e"'^^") for 
small g > and then to extrapolate to g — 1, but typically ($) collapses 
even for small positive g. However, it is evident from our construction 
that Hg is characterized by $o- = 1 for any g < 0, since all the "bad" 
Va{> 0) are replaced by "good" gVa < 0. We can therefore calculate 
{il)g for any g < by a Monte Carlo sampling that is free of the sign 
problem. If {^Tjg is a smooth function of it should then be possible to 
extrapolate to g = 1 (i.e., to the original hamiltonian) from g < 0. We 
emphasize that 5 = is not expected to be a singular point of (il)^; it 
is special o nly in the Monte Carlo evaluation. 

In Fig. ^.3| we exemplify the ^-extrapolation procedure for several 
observables calculated for ^''Fe. In all cases we use polynomial extrapo- 
lations from negative values to the physical case, g — I. The degree of 
the polynomial is usually chosen to be the smallest that yields a per 
degree of freedom less than 1. However, in several studies, like the one 
of the p/-shell nuclei reported in Section HI, we have conservatively cho- 
sen second-order polynomials for all extrapolations, although in many 
cases a first-order polynomial already resulted in x^-values less than 1. 
At T = the variational principle requires that (H) has a minimum at 
g — I. We have incorporated this fact in our extrapolations of ground 
state energies by using a second-order polynomial with zero-derivative 
at 5 = 1. 




(24) 
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2.4 Computational considerations 

The present SMMC code is roughly the fourth major revision of a pro- 
gram whose development began in late 1990 with a single-species, single- 
j-shell version. It is now a modular package of some 10,000 commented 
lines of FORTRAN. All floating point computations are double preci- 
sion (64-bit). 

The package performs all of the functions necessary for shell model 
Monte Carlo calculations: initialization, thermalization of the Metropo- 
lis walk, generation of the Monte Carlo samples, evaluation of static 
observables and response functions (canonical or grand-canonical). Max- 
imum Entropy extraction of strength functions [l^ , and the extrapola- 
tion in g required to solve the sign problem, as discussed in section 2.3. 
The data input and results output use standard shell model conventions, 
so that it is easy to change the two-body matrix elements of the inter- 
action to incorporate additional one- or two-body observables in the 
analysis, or to add or change the orbitals in the calculation. The code 
has been de-bugged and tested extensively against direct diagonalization 
results in the sd- and lower ^/-shells. Its operation by an experienced 
user can be described as "routine," although it takes several weeks to 
acquire that experience. 

Shell model Monte Carlo calculations are extraordinarily well-suited 
to Multiple Inst ruction/ Multiple Data (MIMD) architectures. Very few 
problems are encountered in porting to a new machine, and the opera- 
tion generally takes less than a day. Indeed, our code is "embarrassingly 
parallel" : separate Metropolis random walks are started on each compu- 
tational node, which then produces a specified number of Monte Carlo 
samples at regular intervals during the walk. Data from all of the nodes 
are sent to a central node for evaluation of the Monte Carlo averages 
and their uncertainties. 

To date, we have implemented the parallel version of our code on 
the Intel DELTA and PARAGON machines at Caltech and ORNL (each 
with 512 i860 processors), on the 128-processor IBM SP-1 at ANL, the 
512-processor IBM SP-2 at Maui, and on a Fujitsu VPP500 shared- 
memory vector processor (24 CPU's). In all cases, the ratio of com- 
munications to computation is very low, with efficiencies always greater 
than 95%. 

Table 1 shows benchmarks of our code on various single processors. 
The test calculation involved a canonical ensemble in the full p/-shell 
{Ns = 20 for each type of nucleon, implying 20 x 20 matrices) using a 
realistic interaction. Nt = 32 time slices were used at /3 = 2 MeV~^ 
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Processor 


Peak MF 


Average MF 


Samples/hr. 


i860 


35 


9 


44 


IBM-SP2 Thin66 


56 


36 


179 


ALPHA-400 


56 


28 


141 



Table 1: Benchmarks of the SMMC code in various processors. 



Op 


6 


Is-Od 


12 


Ip-Of 


20 


lp-O.f-Ogg/2 


30 


2s-ld-0g 


30 


2p-l/-0.y9/2-0/i;v2 


44 



Table 2: Matrix dimension for various model spaces. 

(A/3 = 0.0625 McV"^). Thirty static obscrvablcs and seven dynami- 
cal response functions were calculated at a single g-value (see section 
7). Note that the computational speed is independent of the interaction 
and of the number of nucleons occupying the shell. Beyond using library 
subroutines (BLAS and LAPACK routines), no attempt was made to op- 
timize the assembly level code in any of these cases. Approximately 40% 
of the computational effort is in matrix- multiply operations. A signifi- 
cant remainder of the effort goes into building the one-body hamiltonian 
(13%), setting up the one-body evolution operators (15%), and calcu- 
lating two-body observablcs (15%), none of which is easily vectorizable. 
In general, the computation time scales as N^Nf, and is spent roughly 
equally on the dynamical response functions and the static observable 
sampling. 

The memory required for our calculations scales as N^Nf. is 
the average of the squares of the numbers of neutron and proton single- 
particle states. Sample values of Ng for one isospin type are shown for 
various model spaces in Table 2. 

The code is currently structured so that a calculation with Ng = 32, 
six j-orbitals, and Nt = 64 time slices will fit in 12 MB of memory. 
A calculation in the {lp-0f)-{2s-ld-0g) basis has Ng = 50 and would 
require about 64 MB of memory for Nt = 64 time slices and about 
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128 MB for Nt = 128. 



3 Results 

In the past four years SMMC techniques have been applied in a variety 
of ways to nuclei in the p/-shell and other regions using realistic or 
semi-realistic interactions. In this section we will detail several of the 
important results of these calculations. These are, broadly, ground state 
and thermal properties of iron region nuclei, nuclear pair correlations, 
/3/3-decay, and 7-soft studies. 



3.1 Ground state properties of medium mass nuclei 
3.1.1 p/-shell nuclei 

While complete OTiuj calculations can be carried out by direct diagonal- 
ization in the p- and sd-shells, the exponentially increasing number of 
configurations limits such studies in the next (pf) shell to only the very 



lightest nuclei ||25|, 12 . SMMC techniques allow calculation of ground- 
state observables in the full Otko model space for nuclei throughout the 
p/-shell. Here, we discuss a set of such calculations that uses the mod- 
ified KB3 interaction |9|; the single-particle basis is such that Ng = 20 
for both protons and neutrons. These studies were performed for 28 
even-even Ti, Cr, Fe, Ni, and Zn isotopes and 4 odd-odd N = Z nu- 
clei. A more detailed description of the calculations and their results 
are found in Ref . |^ . 

Figure |3.1.l| shows systematic results for the mass defects obtained 
directly from (H). The SMMC results have been corrected for the 
coulomb energy, which is not included in the KB3 interaction, using 

iii 

i/coui = ^'•^ ~ • 0.35 - TTMm + 7T- , 7.289. (25) 

where tt and v are the numbers of valence protons and neutrons, re- 
spectively, and the energy is in MeV. As in Ref. , we have increased 
the calculated energy expectation values by 0.014 • n(n — 1) MeV (where 
n = 7r-|-^ is the number of valence nucleons) to correct for a "tiny" resid- 
ual monopole defect in the KB3 interaction. In general, there is excellent 
agreement; the average error for the nuclei shown is -1-0.45 MeV, which 
agrees roughly with the internal excitation energy of a few hundred keV 
expected in our finite-temperature calculation. 
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FIG 3.1.1 Upper panel (a): Comparison of the mass excesses AM as 
calculated within the SMMC approach with data. Lower panel (b): 
Discrepancy between the SMMC results for the mass excesses and the 
data, SAM. The solid line shows the average discrepancy, 450 keV, 
while the dashed lines show the rms variation about this value (from 



Shell Model Monte Carlo 



17 



Figure 3.1.1 shows the calculated total E2 strengths for selected pf- 



shell nuclei. This quantity is defined as 

B{E2)^{{epQp + e„Qn)^) , (26) 

with 



Y2{9,,cl),), (27) 



where the sum runs over all valence protons (neutrons). The effective 
charges were chosen to be Cp = 1.35 and e„ = 0.35, while we used 
b = 1.01^^/^ fm for the oscillator length. Shown for comparison are the 
B{E2) values for the 0^ — + 2^ transition in each nucleus; in even-even 
nuclei some 20-30% of the strength comes from higher transitions. The 
overall trend is well reproduced. For the nickel isotopes ^8'60^64]\jj^ -j-j^g 
total B{E2) strength is known from (e,e') data and agrees very nicely 
with our SMMC results. 



3.1.2 The 0/5/2lp0gg/2 model space 

In order to investigate heavier systems, the 039/2 was included in the 
pf model space However, we found that the coupling between the 
0/7 /2 and 0^9 /2 orbitals causes significant center-of-mass contamination 
to the ground state, and we therefore close the O/7/2. The model space is 
thus 0/5/2 lp0(79/2. This appears to be a good approximation in systems 
where TV and Z are greater than 28. The monopole terms of this new 
interaction were modified |^ to give a good description of the spectra 
of nuclei in the Ni isotopes. Since ^^Ni is the core of this model space, 
the single-particle energies were determined from the ^^Ni spectrum. 

Shell-model masses must be corrected for coulomb effects and for 
'grand-monopole' terms that are not taken into account when de- 
termining interaction matrix elements from nuclear spectra. Following 
Rcf. ||l^ we have fitted a correction to the calculated binding energies 
which takes into account the coulomb energy and residual effects of the 
monopole terms. The form of this correction is 

He = aiT{TT - 1) -I- Pttv -f 77r -I- 6n{n - 1) -t- en (28) 

where tt, v are the number of valence protons and neutrons, and n = 
TT + ly, and all parameters are in units of MeV. 7 is held fixed at 7.289. 
Minimizing for the difference between experimental and theoretical 
binding energies relative to the ^^Ni core gives a = 0.234, (3 = 0.0156, 
S = 0.0316, and e = 1.828. 



18 



Koonin, Dean, & Langanke 



1250 



1000 



750 



m 500 



250 



48 



. SMMC (total) 
' Expt (0* -> 2*) 
1 Expt (total) 



50 



52 



54 



56 58 
A 



60 



62 



64 



66 



FIG |3.1.1| Comparison of the experimental B{E2,0i ~^ 2f) strengths 
with the total B{E2) strength calculated in the SMMC approach for 
various p/-shell nuclei having either a proton or neutron number of 28. 
For the nickel isotopes 58,60,64]\jj ^^le total B{E2) strength (full squares) 
is known from inelastic electron scattering data (from [p6|). 



Shell Model Monte Carlo 



19 




250 



^ 1.0 - 



N 0.0 - » . Ge * T 

5, ■ -Se ^ . 1 _ 

m ♦Kr - . . MC error 

•c -1 .0 - . A Sr • ,, 
* ^ N=Z 

-2.0 

60 62 64 66 68 70 72 74 76 78 80 82 



FIG 3.1.2 . Calculated mass excesses for nuclei in the mass range A = 
64-80 are compared with experiment (top panel), and the difference 
between experiment and theory is shown (bottom panel). 



Shown in Fig. 3.1.2 



are the experimental and theoretical binding 
energies BE{N,Z) for the nuclei studied here (top panel), and the dif- 
ference between experiment and theory (bottom panel). The overall 
agreement is reasonable. 



3.1.3 Gamow- Teller strengths and distributions 

The Gamow- Teller (GT) properties of nuclei in this region of the peri- 
odic table are crucial for supernova physics The core of a massive 
star at the end of hydrostatic burning is stabilized by electron degener- 
acy pressure as long as its mass does not exceed the appropriate Chan- 
drasekhar mass Mch- If the core mass exceeds Mch, electrons are 
captured by nuclei . For many of the nuclei that determine the elec- 
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tron capture rate in this early stage of the presupernova |3C[ |, Gamow- 
Teller (GT) transitions contribute significantly to the electron capture 
rate. Due to insufficient experimental information, the GT^ transi- 
tion rates have so far been treated only qualitatively in presupernova 
collapse simulations, assuming the GT_|_ strength to reside in a single 
resonance whose energy relative to the daughter ground state has been 
parametrized phenomenologically ; the total GT+ strength has been 



taken from the single particle model. Recent (n,p) experiments |32 -]3q], 
however, show that the GT+ strength is fragmented over many states, 
while the total strength is significantly quenched compared to the single 
particle model. (A recent update of the GT-|_ rates for use in supernova 
simulations assumed a constant quenching factor of 2 @.) 
The total GT strengths are defined as 

S(Gr±) = ((?r±f). (29) 

From Ohuj shell model studies of the GT strengths for sd shell nuclei 
and for light pf shell nuclei it has been deduced that the spin operator 
in (31) should be replaced by an effective one, (Toff = f7/1.26 [T^ . 
This renormalization is not well understood. It is believed to be related 
either to a second-order core polarization caused by the tensor force 
1 37 1, or to the screening of the Gamow- Teller operator by A-hole pairs 
|38|. Using the effective spin operator, our calculated i?(GT4.) are in 
exc ellent agreement with the data deduced from (n,p) reactions (Fig- 
ure 3.1.3| ). Thus, our calculations support both the statement that the 
spin operator is renormalized by a universal factor (1/1.26) in nuclei, 
and the statement that complete shell-model calculations can account 
for the GT strength observed experimentally. 

In a series of truncated shell-model calculations, Aufderheide and 
collaborators have demonstrated that a strong phase space dependence 
makes the Gamow- Teller contributions to the presupernova electron cap- 
ture rates more sensitive to the strength distribution in the daughter 
nucleus than to the total strength |^ . In this work it also became ap- 
parent that complete Ohuj studies of the GT+ strength distribution are 
desirable. Such studies are now possible using the SMMC approach. 

To determine the GT4. strength distribution, we have calculated the 
response function of the ut^ operator, Rgt(j), as defined in Eq. (|^). 
As the strength function Sgt{E) is the inverse Laplace transform of 
Rqt{t), we have used the Maximum Entropy technique, described in 
Ref. |l8|, to extract S'gt(£^)- 

As first examples we have studied several nuclei (^^V, ^^'^^Fe, ^^Mn, 
58,60,62,64 j^j^ and ^^Co), for which the Gamow-Teller strength distribu- 
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FIG 3.1.3. Comparison of the renormalized total Gamow- Teller 
strength, as calculated within the present SMMC approach, and the 
experimental B{GT+) values deduced from {n,p) data [|2|-|3^. [For 
the odd nuclei (^^V, ^^Mn and ^^Co) the SMMC calculations have been 
performed at /3 = 1 MeV^^ to avoid the odd-j4 sign problem. For even- 
even nuclei the GT strength calculated at this temperature is still a good 
approximation for the ground state value, as discussed in the text.] 
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tion in the daughter nucleus is known from {n,p) experiments |32 -|3 



However, note that electron capture by these nuclei plays only a minor 
role in the presupernova collapse. As SMMC calculates the strength 
function within the parent nucleus, the results have been shifted using 
the experimental Q-values, and the Coulomb correction has been per- 
formed using Eq. (^5|). For all nuclei, the SMMC approach calculates 
the centroid and width of the strength distribution in good agreement 
with data (see Fig. 3.1.3). The centroid of the GT+ strength distribu- 
tions is found to be nearly independent of temperature, while its width 
increases with temperature. 

Following the formalism described in Refs. ||3l|, the Gamow- 
Teller contributions to the electron capture rates under typical presu- 
pernova conditions have been calculated, assuming that the electrons 
have a Fermi-Dirac distribution with a chemical potential adopted from 
the stellar trajectory at the electron-to-nucleon ratio corresponding to 
the particular nucleus |p(| . The calculations have been performed using 



both the SMMC and experimental GT_|_ strength distributions 1 32 



The two electron capture rates agree within a factor of two for tempera- 
tures T — (3 — 5) X 10^ Kelvin, which is the relevant temperature regime 



in the presupernova collapse |30 . Thus, for the first time it is possible to 
calculate, with reasonable accuracy, the electron capture rate for nuclei 
like ^^Co or ^^Ni, which dominate the electron capture process in the 



early presupernova collapse 1 30 



3.2 Pair correlations 

The residual nuclear interaction builds up pairing correlations in a nu- 
cleus. Introducing nucleon creation operators a^, these correlations can 
be studied by defining pair creation operators 



1 



\/l + Sab L JM 

for proton-proton or neutron-neutron pairs, and 



(30) 



1 



± 



\/2(l + 5ab) 

'1 }• 



JM 



(31) 
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. FIG |3.1.3[ Shown are SMMC GT+ strength distributions (solid line) 
for various nuclei in the iron region. The energies refer to the daugh- 
ter nuclei. The dashed histograms show the experimental strength 
distribution as extracted from {n,p) data. (Following Q the calcu- 
lated strength distributions have been folded with Gaussians of width 
1.77 MeV to account for the experimental resolution.) 
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for proton-neutron pairs where "+(—)" is for T = 0(T = 1) pn-pairing. 
With these definitions, we construct a pair matrix 

ML, = J2 (Am Ua , Jb)AjM (jc, Jd)) , (32) 

M 

where a = {ja,jb} and a' = {jc,jd} and the expectation value is in the 
ground state or canonical ensemble at a prescribed temperature. The 
pairing strength for a given J is then given by 

P(J)= ^ Afi,„,. (33) 

a>a' 

An alternative measure of the overall pair correlations is given in 
terms of the BCS pair operator 

Ak=E4M(«)- (34) 

a 

The quantity ^j(^(A j^^Ajj\/) is then a measure of the number of nu- 
cleon pairs with spin J. For the results discussed below, the BCS- 
like definition for the overall pairing strength yields the same qual- 
itative results for the pairing content as the definition (p3|). Some 
SMMC results for BCS pairing in nuclei A = 48 — 60 are published 
in Refs. |6[|o[|4ll ||. 



With our definition in (|3^) the pairing strength is positive at the 
mean-field level. The mean-field pairing strength, PupiJ), can be de- 
fined by replacing the expectation values of the two-body matrix ele- 
ments in the definition of M"^ by 

(a|4«3a4) riin2 {S13S24 - (523^14) , (35) 

where Uk — (ajjOfc) is the occupation number of the orbital k. This 
mean-field value provides a baseline against which true pair correlations 
can be judged, by defining 

-Pcor(J) = P(J)-Pmf(J) . (36) 



3.2.1 Ground state pair correlations 

In a first project we studied the pair correlations in several even-even 
pf shell nuclei (54>56,58pg g^^^ ^^Cr) As expected, we found a large 
excess of J = 0+ like-particle pairing in the ground states of these 
nuclei. With increasing temperature, these pairing correlations decrease 
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and at around T = 1 MeV the like-particle pairs break in these nuclei 
(pairing phase transition). Our calculations also indicate that isoscalar 
proton-neutron (mainly J — 1+) pairs persist to higher temperatures. In 
Ref. we have related the thermal dependence of several observables 
to the temperature dependence of associated pairing correlations. 

It has long been anticipated that J = 0+ proton- neutron correlations 
play an important role in the ground states oi N = Z nuclei. These 
correlations were explored with SMMC ior N ^ Z nuclei with A = 
48 - 58 in the p/-shell and ^4 64 - 74 in the 0/5/2lpOgg/2 space. 
As the even-even N ^ Z nuclei have isospin T — 0, {A^A) is identical in 
all three isovector 0+ pairing channels. This symmetry does not hold for 
the odd-odd N — Z nuclei in this mass region, which usually have T = 1 
ground states, and {A''' A) can differ for proton-neutron and like-nucleon 
pairs. (The expectation values for proton pairs and neutron pairs are 
identical.) 

We find the proton-neutron pairing strength is significantly larger 
for odd-odd N = Z nuclei than in even-even nuclei, while the 0^ proton 
and neutron pairing shows the opposite behavior, in both cases leading 
to a noticeable odd-even staggering, as displayed in Fig. 3.2.1 for the 
pf shell. Due to the strong pairing in the /7/2 orbital, all three isovec- 
tor 0^ channels of the pairing matrix exhibit essentially only one large 
eigenva lue w hich is used as a convenient measure of the pairing strength 
in Fig. 3.2.1 . This staggering is caused by a constructive (destructive) 
interference of the isotensor and isoscalar parts of A^A in the odd-odd 
(even-even) N = Z nuclei. The isoscalar part is related to the pairing 
energy and is roughly constant for the nuclei we have studied. 

Fig. [3.2.1 shows the correlated pairs for N = Z nuclei in the A = 64— 
74 region of the 0/5/2lp0(7g/2 space. The correlated pairs exhibit a strong 
J = 0, T — 1 like-particle staggering for the even-even and odd-odd 
N — Z systems, while the number of correlated proton-neutron pairs 
is much larger than the like-particle number for the odd-odd systems. 
The correlated pairing behavior ofN = Z,N:^Z + 2, N = Z + 4: 
nuclei is shown in Fig. 3.2.1, where one clearly sees the decrease in T=l 
proton-neutron pairing as one moves away from N = Z. As discussed 
by Engel et al. , increasing the number of neutron pairs increases the 
collectivity of the neutron condensate, making fewer neutrons available 
to pair with protons. As a result, the protons pair more often with one 
another, although their number has not changed, and the np pairing 
drops drastically. A weak neutron shell closure is evident in the Kr 
isotopes as the neutrons fill to N=40 at ^^Kr. Although the neutron 
pairing correlations decrease here, they are not zero as the 179/2 has 
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Fig 3.2.1. Largest eigenvalues for the J = 0, T = 1 proton-proton 
(left) and proton-neutron (right) pairing matrix as a function of mass 
number. 
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Fig 3.2.1 . Correlated J = 0, T ~ 1 pairs in the proton- proton and 
proton- neutron channels for selected isotope chains in the 0/5/2lp0(7g/2 
model space. 



some occupation. 



3.2.2 Pairing and rotation 

A recent study [^of the 7 decays of the odd-odd N — Z nucleus ^^Rb 
revealed the phenomenon of an isospin band crossing at modest excita- 
tion energies. While the ground state rotational band can be identified 
as being formed from the T — 1 isobaric analogue states of ^''Kr, a 
T = band becomes energetically favored over the T = 1 band with 
increasing rotational frequency. To study this isospin band crossing, 
we have performed a cranked SMMC calculation in which the shell 
model hamiltonian is replaced by — > -f ojJz- Note that since Jz is 
a time-odd operator, the sign problem is reintroduced, and good statis- 
tical sampling imposes a limit on uj. For the calculations presented in 
Table 3, the largest cranking frequency was w = 0.4 MeV. In agreement 
with experiment, we find T — 1 for the ground state (w = 0). How- 
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ever, with increasing frequency T = states are mixed in (and even 
dominate at large angular momenta) and (T^) decreases, as can be seen 
in Table 3. Thus our calculation confirms that ^^Rb changes from a 
T = 1 dominated system to one dominated by T = states with in- 
creasing rotational frequency. To understand the apparent crossing of 
the T = 1 and T = bands, we have studied the various pair correla- 
tions as a function of rotational frequency. We find that the isovector 
J = correlations and the aligned isoscalar J = 9 pn correlations are 
most important in this transition. These correlations are plotted in 
Fig. 3.2.2 as a function of (Jz), where we have defined pair correlations 
by Eq. (t36[). Strikingly, the largest pair correlations are found in the 
isovector J = and isoscalar J — 9 proton-neutron channels at low 
and high frequencies, respectively. Furthermore, the variation of isospin 
with increasing frequency reflects the relative strengths of these two pn 
correlations. Our calculation clearly confirms that proton-neutron cor- 
relations determine the behavior of the odd-odd N = Z nucleus ^"^Rb, 
as already supposed in Q. 

With increasing frequency, the isovector J = pn correlations de- 
crease rapidly to a constant at (Jz) « 3. This behavior is accompanied 
by an increase of the pn correlations in the maximally aligned channel, 
J = 9, T = which dominates ^''Rb at rotational frequencies where the 
J — pn correlations become small. Furthermore, although J = 8, 
T = 1 pairs exist, th ey do not exhibit correlations beyond the mean 
field, as shown in Fig. 3.2.2 , We note that the experiment also indicates 
that the mixing of T = states sets in near J = 3. 



3.2.3 Pairing and temperature 

Another striking difference in proton-neutron pairing can be found in 
the thermal properties of odd-odd (e.g. ^°Mn) and even-even (e.g. ^^Fe) 
N = Z nuclei |4l[. As in other even-even N = Z nuclei, the ground 
state isospin of ^^Fe is T = and isospin symmetry forces pp, nn, 
and np pairing to be identical. With increasing temperature, T = 1 
components are slowly mixed in, breaking the symmetry between nn-pp 
and np pairing. However, (T^) « for T < 1 MeV and the symmetry 
holds. As J = 0^ pairs break at around this temperature in even-even 
nuclei in this mass range, the phase transition is clearly noticeable in all 
three isovector pairing channels. 

As discussed above, the ground state of the odd-odd nucleus ^°Mn 
(with isospin T = 1) is dominated by 0"*" pn correlations. As a striking 



feature shown in Fig. 3.2.3, the proton-neutron pairing decreases rapidly 
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<J,> 



Fig 3.2.2. Selected pair correlations and the proton occupation 
number (bottom panel) as a function of {Jz)- The top panel shows the 
isovector J = pp and pn correlations, while the isoscalar J = 9 and 
isovector J = 8 pn correlations are shown in the middle panel. 
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Table 3: Summary of various properties of ^^Rb. These results, quoted 
as a function of the cranking frequency lo, include (T^), the pro- 

ton-proton J — Q pairing correlations (J = 0, pp), J — 0, T — 1, pn 
pairing correlations (J = 0, pn). 



uj MeV 




(T2) 


J = pp 


J = pn 


0.00 


3.5 ± 0.6 


1.97 ± 0.05 


0.21 ± 0.05 


1.81± 0.08 


0.10 


2.1 ± 0.9 


1.825± 0.03 


0.11 ± 0.05 


1.53±0.08 


0.20 


2.9 ± 0.7 


1.98 ± 0.10 


0.37 ± 0.06 


1.26±0.08 


0.25 


2.6 ± 1.1 


1.72 ± 0.16 


0.45 ± 0.06 


0.60±0.08 


0.30 


3.3 ± 0.8 


1.77 ± 0.13 


0.46 ± 0.06 


0.52±0.08 


0.35 


6.6 ± 3.1 


1.40 ± 0.14 


0.46 ± 0.06 


0.44±0.08 


0.40 


11.3 ± 1.9 


0.69 ± 0.20 


0.48 ± 0.06 


0.48±0.08 



with temperature and has dropped to the mean-field value by T ~ 
0.75 MeV, while the like-particle pairing remains roughly constant to 
1.1 MeV. The vanishing of the pn correlations is accompanied by a 
change in isospin, which decreases from the ground state value (T^) = 2 
to (T^) w 0.2 at temperatures near T = 1 MeV. 

As required by general thermodynamic principles, the internal en- 
ergy increases steadily with temperature. The heat capacity C(T) = 
dE/dT, with E — (H), is usually associated with the level density pa- 
rameter a by C{T) — 2a{T)T. As is typical for even-even nuclei |^ 
a(T) increases from a = at T=0 to a roughly constant value at tem- 
peratures above the phase transition. We find a{T) « 5.3 ± 1.2 MeV~^ 
at T > 1 MeV, in agreement with the empirical value of 6.5 MeV^"'^ |Q 
for ^^Fe. At higher temperatures, a(T) must decrease due to the finite 
model space of our calculation. 

The temperature dependence of E in ^°Mn is significantly different 
from that in even-even nuclei. As can be seen in Fig. 3.2.3| , E increases 
approximately linearly with temperature, corresponding to a constant 
heat capacity C{T) w 5.4 ± 1 MeV"-^; the level density parameter de- 
creases like a(r) ~ in the temperature interval between 0.4 MeV 
and 1.5 MeV. We note that the same linear increase of the energy with 
temperature is observed in SMMC studies of odd-odd N — Z nuclei 
performed with a pairing-|-quadrupole hamiltonian | |4^ . 

To investigate the apparent differences in the thermal behavior of 
even-even and odd-odd N — Z nuclei in a more systematic way, we 



Shell Model Monte Carlo 



31 




1.0 
T (MeV) 




1.0 2.0 

T (MeV) 



Fig. 3.2.3. Thermal properties of ^ Fe. The SMMC results are shown 



with error bars, while the lines indicate the mean-field values for the 
respective pair correlations. 
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have performed SMMC calculatfons for ^''Mn and 5°Cr using a pair- 
ing+quadrupole hamiltonian. Note that the quadrupole part aUows for 
T = pn correlations. The strength of the quadrupole component 
has been adjusted to reasonably reproduce SMMC results for the re- 
alistic KB3 interaction. This hamiltonian has no sign problem, thus 
significantly reducing the statistical uncertainties (and also potential 
systematic errors related to the g-extrapolation). Despite its simplicity, 
the hamiltonian nevertheless embraces much of the essential degrees of 
freedom governing the thermal response at low temperatures. 

The ground states of ^"Cr and ^°Mn belong to the same T = 1 mul- 
tiplet and are therefore isobaric analogues. This is correctly recovered 
in the SMMC calculation as the energy expectation values at low tem- 
peratures are, within error bars, identical. However, (T^) = 1.89 ±0.07 
for ^°Mn at r = 0.33 MeV (the lowest temperature at which we have 
performed these SMMC calculations) indicates some isoscalar compo- 
nents being mixed in. (Experimentally, the lowest T = state in 
^''Mn is at an excitation energy of 0.2 MeV.) With increasing tempera- 
ture the relative strength of the isovector component in ^"Mn weakens. 
Correspondingly, (T^) decreases to 1.34 ± 0.04 at T = 0.8 MeV. As 
isospin T = states cannot be formed in ^°Cr due to the neutron ex- 
cess ((T^) — 2.1 at T = 0.8 MeV), this nucleus has fewer degrees of 
freedom than ^"Mn. As a consequence ^°Mn has a higher excitation 
energy, or relatedly a higher level density, than ^"Cr at modest tem- 
peratures (at T = 0.8 MeV (H) = -7.55 ± 0.13 MeV for ^OMn and 
—8.28 ± 0.08 MeV for ^°Cr). At ev en hig her temperatures, the like- 
particle correlations break (see Fig. 3.2.3 ). For both nuclei, isospin 
states with T > 2 get mixed in and (T^) increases. For temperatures 
above the phase transition (T > 1 MeV), the thermal properties for 
both nuclei {H, ,Ql,Qp,Qn^ the latter two with the appropriate 

scaling by the proton and neutron numbers) become almost the same. 
At T < 1 MeV, the isovector 0+ pairing correlations behave as for the 
realistic KB3 interaction. In particular, there is the noticeable excess 
of pp and nn correlations in the even-even nucleus ^"Cr at low temper- 
atures, while the odd-odd ^'^Mn is dominated by pn correlations. The 
latter decrease rather rapidly with temperature, while the pp and nn 
correlations in ^"Cr show the behavior characteristic of a phase transi- 
tion near 1 MeV, which is typical for even-even nuclei in this mass range 
studied with the KB3 interaction. 

However, there are a few differences between calculations with the 
simple hamiltonian and the realistic one. First, for the simple hamil- 
tonian all isovector correlations decrease much more slowly at temper- 
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atures exceeding 1 MeV. The origin for this is the missing isoscalar pn 
correlations, which dominate nuclear properties at T > 1 MeV |4^. 
Their persistence at these higher temperatures further suppresses the 
isovector correlations when the realistic hamiltonian is used. Second, 
quantities that are sensitive to np correlations are not correctly de- 
scribed by the simple hamiltonian. Among these are the Gamow- Teller 
strengths. Here the realistic hamiltonian yields B{GT+) = 5.2 ± 1.8 for 
^"Mn and 2.2±0.2 for ^°Cr (after the renormalization of the spin opera- 
tor). In contrast, for the simple hamiltonian we find B(GT+) = 4.5±0.1 
for ^°Mn and 3.8±0.01 for ^"Cr. These values have to be compared with 
the independent-particle model estimates: 6.9 (^°Mn) and 5.1 (^°Cr). 
While the large error bars make a comparison for ^'^Mn meaningless, the 
isovector pn pairing accounts for somewhat less than half of the GT+ 
quenching in ^''Cr, the reminder being due to isoscalar pn correlations 
between spin-orbit partners, introduced by the (ctt)^ piece missing in 
the simple hamiltonian. Third, the level density at moderate temper- 
atures (T « 1 MeV) is less for the realistic hamiltonian than for the 
simple one; one of the reasons again is the missing isoscalar pn correla- 
tions that push levels to higher energies. 

Nevertheless these SMMC calculations with the sign-problem-free 
hamiltonian are quite illustrative and can be extended to odA-A and 
odd-odd N ^ Z nuclei. First calculations show that one can investigate 
these nuclei down to temperatures of order 0.5 MeV, before the sign- 
problem seriously enters again. 



The second-order weak process [Z, A) ^ {Z + 2, A) + 2e^ + 2De is an 
important "background" to searches for the lepton-number violating 
neutrinoless mode, {Z,A) — > (Z -I- 2, A). The calculation of the nuclear 
matrix element for these two processes is a challenging problem in nu- 
clear structure, and has been done in a full pf model space for only the 
lightest of several candidates, ''^Ca. P.B. Radha et al. have performed 
first Monte Carlo calculations of the 2^ (3f3 matrix elements in very large 
model spaces Q. 

In two-neutrino double /3-decay, the nuclear matrix element of inter- 
est is 



where |io) and |/o) are the 0+ ground states of the initial and final even- 



3.3 /5/3-decay 




(37) 
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even nuclei, and |m) is a l''^ state of the intermediate odd-odd nucleus; 
the sum is over all such states. In this expression, G = (TT_ is the 
Gamow- Teller operator for /3~-decay (i.e., that which changes a neutron 
into a proton) and f2 = (Ei^ + Ef^)/2. A common approximation to 
M'^'^ is the closure value, 

M^'' = ^ (38) 

where E is an average energy dominator and 

Me = ^(/o|G|m)(m|GKo) = (/o|G • G|zo) • (39) 

in 

SMMC methods can be used to calculate both Mc and M"^" . To do 
so, consider the function 



(T,r') = (e-^(^+^')G+ • Gte--f^^Ge-^^'G) 



Z 



-(/3-r-r')ffGt . Gte-«Ge--'^G 



, (40) 



where Z = Tr^ e^^^ is the partition function for the initial nucleus, H is 
the many-body hamiltonian, and the trace is over all states of the initial 
nucleus. The quantities (/3 — r — r') and r play the role of the inverse 
temperature in the parent and daughter nucleus respectively. A spectral 
expansion of shows that large values of these parameters guarantee 
cooling to the parent and daughter ground states. In these limits, we 
note that (/i(r, r' = 0) approaches e~'^'5|Mcp, where Q — Ef ~ E'j is 
the energy release, so that a calculation of 0(t, 0) leads directly to the 
closure matrix element. If we then define 

7?(T,r)^ f dT'cj,[T,T')e-^'QI\ (41) 
Jo 

and 

"^'<^-)-^' (*^' 

it is easy to see that in the limit of large t, (/3— r — t'), and T, M^''(T, r) 
becomes independent of these parameters and is equal to the matrix 
element in Eq. (^. 

In the first applications, Radha et al. calculated the 2v matrix ele- 
ments for "^^Ca and ^^Ge [Q. The first nucleus allowed a benchmarking 
of the SMMC method against direct diagonalization. A large-basis shell 
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model calculation for ''^Ge has long been waited for, as '^Ge is one of 
the few nuclei where the 2i>P(3 decay has been measured precisely and 
the best limits on the decay mode have been established ^ |5^ . 

To monitor the possible uncertainty related to the ^-extrapolation 
in the calculation of the 2i/ matrix element for ^^Ge, SMMC studies 
have been performed for two quite different families of sign-problem- 
free hamiltonians {x = oo and X = 4). The calculation comprises the 
complete (O/5/2, Ip, Ogg/2) niodel space, which is significantly larger than 
in previous shell model studies [^. The adopted effective interaction 
is based on the Paris potential and has been constructed for this model 
space using the Q-box method developed by Kuo |Q . 

As is shown in Fig. upon linear extrapolation, both families of 
hamiltonians predict a consistent value for the 21/ matrix element of 
^^Ge. The results M^'' = 0.12 ± 0.07 and M^" = 0.12 ± 0.06 are only 
shghtly lower than the experimental values (M^" = 0.22 ± 0.01 
This comparison, however, should not be overinterpreted, as the detailed 
reliability of the effective interaction is still to be checked. 

It is interesting that the closure matrix element found in the SMMC 
calculation and the average energy denominator {Mc — —0.36 ± 0.37, 
E = -3.0 ± 3.3 MeV and = 0.08 ± 0.17, E = 0.57 ± 1.26 MeV for 
the two families of hamiltonians with x = 00 and x = 4, respectively) 
are both significantly smaller than had been assumed previously. This 
is confirmed by a recent truncated diagonalization study 



3.4 7-soft nuclei 

Nuclei with mass number 100 < A < 140 are believed to have large 
shape fluctuations in their ground states. Associated with this softness 
are spectra with an approximate 0(5) symmetry and bands with energy 
spacings intermediate between rotational and vibrational. In the geo- 
metrical model these nuclei are described by potential energy surfaces 
with a minimum at /? ^ but independent of 7 . Some of these nu- 
clei have been described in terms of a quartic five-dimensional oscillator 
iQ. In the Interacting Boson Model (IBM) they are described by an 
0(6) dynamical symmetry |57| |58| , [59[| . In the following we review the 
first fully microscopic calculations for soft nuclei with 100 < A < 140 




For the two-body interaction we used a monopole (J = 0) plus 
quadrupole (J = 2) force ^l| supplemented by a collective quadrupole 
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The 2^ matrix element for ^^Ge calculated within SMMC 
studies based on two families of hamiltonians which are free of sign 
problems. The physical values are obtained by linear extrapolation to 
g = 1. The experimental value for this matrix element is indicated 
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interaction: 

Ho = 



E - h ■■ T.(-rQ,Q-, ■■ . (43) 



, 2A+1 



where :: denotes normal ordering. The single particle energies and the 
other parameters were determined as described in Ref. [ |60| . 

We begin discussion of our results with the probability distribu- 
tion of the quadrupole moment. This is obtained from the shape of 
each Monte Carlo sample, including the quantum-mechanical fluctua- 
tions through the variance of the Q operator for each sample, = 
Tr(Q2[/_^)/Tr([/^) - (Q)^. The shape distribution P(/3,7) can be con- 
verted to a free energy surface as discussed in Ref. [|o| . 



The shape distributions of ^^^Te and ^^''Xe are shown in Fig. 3.4 at 
different temperatures. These nuclei are clearly 7-soft, with energy min- 
ima at /3 ~ 0.06 and f3 ^ 0.15, respectively. Energy surfaces calculated 
with Strutinsky-BCS using a deformed Woods-Saxon potential also 
indicate 7-softness with values of /? comparable to the SMMC values. 
These calculations predict for -'^^^Xe a prolate minimum with /3 ~ 0.20 
which is lower than the spherical configuration by 1.7 MeV but is only 
0.3 MeV below the oblate saddle point, and for ^^^Te a shallow oblate 
minimum with /3 « 0.03. These 7-soft surfaces are similar to those 
obtained in the 0{6) symmetry of the IBM, or more generally when 
the hamiltonian has mixed C/(5) and 0(6) symmetries but a common 
0(5) symmetry. In the Bohr hamiltonian, an 0(5) symmetry occurs 
when the collective potential energy depends only on /3 [^. The same 
results are consistent with a potential energy V{f}) that has a quartic 
anharmonicity p^ , but with a negative quadratic term that leads to a 
minimum at finite /3. 

The total E2 strengths were estimated from (Q^) where Q = CpQp + 
GnQn is the electric quadrupole operator with effective charges of Cp = 
1.5e and e„ = 0.5e, and B{E2; 0^ 2^) determined by assuming that 
most of the strength is in the 2f state. Values for B{E2;0^ — > 2^) of 
663 ± 10, 2106 ± 15, and 5491 ± 36 e^fm'* were found, to be compared 



with the experimental values |63 of 1164, 3910, and 9103 e'^fm^ for 
i24Sn, i28Te, and ^^^Xe, respectively. Thus, the SMMC calculations 
reproduce the correct qualitative trend. The 2f excitation energies were 
also calculated from the E2 response function. The values of 1.12 ± 
0.02, 0.96 ± 0.02, and 0.52 ± 0.01 MeV are in close agreement with the 
experimental values of 1.2, 0.8, and 0.6 MeV for ^^^Te, and ^^-^Xe, 

respectively. 
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FIG 3.4. Contours of the free energy (as described in the text) in the 
polar-coordinate /? — 7 plane for ^^®Te and ^^*Xe. Contours are shown 
at 0.3 MeV intervals, with lighter shades indicating the more probable 
nuclear shapes (from jGOf). 
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Another indication of softness is the response of the nucleus to ro- 
tations, probed by adding a cranking field uJ^ to the hamiltonian and 
examining the moment of inertia as a function of the cranking frequency 
oj. For a soft nucleus one expects a behavior intermediate between a 
deformed nucleus, where the inertia is independent of the cranking fre- 
quency, and the harmonic oscillator, where the inertia becomes singular. 



This is confirmed in Fig. 3.4 which shows the moment of inertia I2 for 
^^"'Xe and ^^^Te as a function oiu, and indicates that ^^*Te has a more 
harmonic character. The moment of inertia for w = in both nuclei is 
significantly lower than the rigid body value (« 43?i^/MeV for A — 124) 
due to pairing correlations. 



Also shown in Fig. 3.4 are {Q ) where Q is the mass quadrupole, 
the BCS-like pairing correlation (A^A) for the protons and {Jz) (neu- 
tron pairing is less affected, and therefore not shown). Notice that the 
increase in I2 as a function of oj is strongly correlated with the rapid 
decrease of pairing correlations and that the peaks in I2 are associated 
with the onset of a decrease in collectivity (as measured by (Q^)). This 
suggests band crossing along the yrast line associated with pair breaking 
and alignment of the quasi-particle spins aX uj k, 0.2 MeV {{Jz) ~ IK) 
for i28Te and w w 0.3 MeV {{Jz) « ll?i) for i24xe. The results are 
consistent with an experimental evidence of band crossing in the yrast 
sequence of ^^"^Xe around spin of 10 fi |6^. The alignment effect is 
clearly seen in the behavior of (Jz) at the lower temperature, which 
shows a rapid increase after an initial moderate change. Deformation 
and pairing decrease also as a function of temperature. 



The total number of J-pairs (nj = riaj) in the various pairing 
channels was also calculated. Since the number of neutrons in ^^"^Xe is 
larger than the mid-shell value, they are treated as holes. For J = 
and J — 2, one can compare the largest Uaj with the number of s and d 
bosons obtained from the 0(6) limit of the IBM. For ^"^^Xe the SMMC 
(IBM) results in the proton-proton pairing channel are 0.85 (1.22) s 
(J = 0) pairs, and 0.76 (0.78) d {J — 2) pairs, while in the neutron- 
neutron channel we find 1.76 (3.67) s pairs and 2.14 (2.33) d pairs. For 
the protons the SMMC c? to s pair ratio 0.89 is close to its 0(6) value 
of 0.64. However, the same ratio for the neutrons, 1.21, is intermediate 
between 0(6) and SU{i) (where its value is 1.64) and is consistent with 
the neutrons filling the middle of the shell. The total numbers of s and 
d pairs — 1.61 proton pairs and 3.8 neutron (hole) pairs — are below the 
IBM values of 2 and 6, respectively. 
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4 Illustrative calculations 

4.1 Giant dipole resonances 

To describe the GDR in ^^O, the Op-lsO(i-0/7/2 shell model space was 
used | |65[ |, and the residual interaction was composed of pairing and 
multipole terms with the specific parametrization given in jBSf , 

^rcs = Xi(9p - Qnf + X2{ql + qlf + X4{qt + qtf ■ (44) 



Shown in Fig. 11 are the results obtained for ^^O for the dipole 
transitions D = l/2{q^-qD at temperatures T =0.5, 1.0, and 2.0 MeV. 
Although these calculations require further refinement, some general 
remarks are warranted. First, from the values of i?(r) at r = 0, the 
total dipole strength is roughly constant as a function of temperature, 
as is the first moment of the strength function. The width, or second 
moment, cannot be reliably extracted from the calculations presented 



in Fig. 4.1 due to large Monte Carlo errors. Still, the calculations are 
representative of what may be obtained with more samples (yielding 
smaller errors). A future study will involve a more complete model 
space description, and implementation of matrix-multiply stabilization 
so that lower temperatures may be studied. 



4.2 Multi-major shell calculations 

Numerous problems in nuclear physics would benefit from the ability to 
calculate shell model observables in more than a single oscillator shell. 
For example, the structure of light neutron-rich nuclei in the sc?-shell 
requires inclusion of /p-shell orbitals in order to properly describe their 
deformation and cross shell characteristics. Experimentally, this region 
will be probed to understand the nature of weakly bound systems, and 
the response of nuclei near the neutron drip line to various physical 
probes. Multi-?ia; SMMC calculations should be no more difficult com- 
putationally than the Ohu calculations that are now routine, although 
they do require increased memory and computational cycles. 

When more than one oscillator shell is included, a new complication 
arises: excitations of the center-of-mass (CM) become mixed into the 
real nuclear excitations. This leads to spurious components of wave func- 
tions being incorporated into the calculations of observables. We have 
investigated two methods of removing these spurious components in the 
SMMC. Projecting out these components leads to another sign problem 
and has proven unworkable. The Gloeckner-Lawson prescription Bq] has 
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FIG 4.1 Response function (left panels) and their excitation strength 
distributions (right panels) for isovector dipole transitions in ^^O at 
T = 0.5, 1.0, and 2.0 MeV. 
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been shown to work for fairly small CM multipliers (< 5) in comparisons 
with direct diagonalization calculations. Unfortunately, this method re- 
quires a "complete" nhui space to be fully effective and SMMC meth- 
ods cannot work in such a complete space The Gloeckner-Lawson 
prescription will still remove the spurious solutions at the cost of also 
slightly affecting the non-spurious solutions. The best value of the mul- 
tiplier is currently being investigated. Values of order 1 seem to be 
best. 

Our current research focuses on the iV = 20 and 28 shell closures for 
very neutron-rich nuclei where we may investigate how exact the shell 
closure really is for nuclei such as ■^^Mg and ''^Ar. (In the sdfp shell this 
calculation allows excitations of up to 16huj for •^^Mg.) Experimentally, 
■^^Mg is highly deformed An sc?-shell calculation gives roughly 30% 
of this deformation, while the inclusion of the p/-shell yields much of 
the remaining deformation due to two-particle excitations into the pf- 
shell. There is no well-established interaction for the sdpf shell, so we 
have chosen to use the WBMB interaction of Ref. |^. In order to 
find agreement with the experimental mass excess the single particle 
energies of the pf shell were lowered by 4 MeV relative to the sd shell. 
We are interested in both the deformation of '^^Mg, and the neutron-rich 



sulfur isotopes which have recently been experimentally measured] 70 



In Fig. 4.2 we show for several nuclei in the region, the calculated mass 
excess for this interaction compared to experiment. We also show the 
total calculated B{E2) vs. the experimentally measured B{E2) to the 
first 2+ state. The final panel shows the expectation value of the center- 
of-mass hamiltonian. For all of the nuclei, a spurious 2huj solution would 
have (Hcm) ~ 20 MeV. The centcr-of-mass contamination is seen to be 
quite small. 



5 Concluding remarks 

We have demonstrated in the preceding two sections the powers and 
limitations of the SMMC technique in describing nuclear properties. 
The SMMC technique, while often non-trivial to apply to a given sit- 
uation, indeed sheds light on various aspects of nuclear systems. The 
most important studies mentioned above include the quenching of the 
Gamow- Teller strength, thermal and rotational properties of pair corre- 
lations, and the evolution of nuclear shapes with temperature. 

We have witnessed enormous computational advances in the last sev- 
eral years. Large-scale parallel computing for scientific problems contin- 
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ues to move forward quickly. Thus some of the more difhcult problems 
that have not been treated by SMMC at present could be pursued in the 
near future. Among these is the GDR in a several fiui space; studies of 
nuclei in the Pb region; studies of weakly bound systems; and studies at 
very low temperatures [requiring matrix-multiply stabilization discussed 
in Ref. Q]. The low temperature studies would allow for a detailed 
description of strength functions not currently available to us. The ad- 
vent of radioactive beam experiments will also allow for interesting new 
physics to be studied within the realm of the shell model. 

The ground state and thermal properties of nuclear matter are an- 
other intriguing application of SMMC methods. One approach is to use 
single particle states that are plane waves with periodic boundary condi- 
tions, and a G-matrix derived from a realistic inter- nucleon interaction; 
the formalism and algorithms we have presented here are then directly 
applicable. An alternative approach is to work on a regular lattice of 
sites in coordinate space and employ Skyrme-like effective interactions 
that couple neighboring sites; the calculation is then similar to that for 
the Hubbard model for which special techniques must be used to handle 
the large, sparse matrices involved ||2l|. Both approaches are currently 
being pursued. 

Otsuka and collaborators have recently proposed a hybrid scheme 
whereby some of the SMMC sampling methods are used to select a 
many-body basis, which is then employed in a conventional diagonal- 
ization. The brunt of their effort involves finding the local Hartree- 
Fock minima of the system, and then sampling around those minima. 
J-projection is approximately incorporated. With this approach one is 
able to obtain the lowest few discrete levels and transitions. The method 
is a way of truncating the many-body basis before diagonalization, and 
is thus complementary to SMMC in the same way as the conventional 
shell model. 

The work discussed above is the culmination of the efforts of many 
people including Y. Alhassid, C. Johnson, G. Lang, E. Ormand, P.B. Radha, 
P. Vogel, and more recently M.T. Ressell, and J. White. This work 
was supported in part by the National Science Foundation, Grant Nos. 
PHY94-12818 and PHY94-20470. Oak Ridge National Laboratory (ORNL) 
is managed by Lockheed Martin Energy Research Corp. for the U.S. De- 
partment of Energy under contract number DE-AC05-96OR22464. DJD 
acknowledges an E.P. Wigner Fellowship from ORNL. Compuational 
resources were provided by the Center for Advanced Computational 
Research at Caltech, the Maui High Performance Computing Center, 
The RIKEN computer center, and Center for Computational Sciences 
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